Verification and prediction of timing behavior of circuit designs is an integral part of electronic design automation (EDA). This is accomplished through computer aided circuit analysis programs, generally referred to as circuit simulators. Simulating a circuit's behavior prior to manufacturing provides insights into the behavior of electronic circuit designs. Accurate prediction of timing behavior during various stages of the design cycle is necessary to avoid expensive post-design corrective measures, delayed time-to-market and failed design cycles.
The circuit is typically described by the connectivity information among various components and devices, using various means, such as netlists, schematic diagrams, textual data etc. The devices/components are represented by appropriate models representing their physical behavior. Most models try to replicate the behavior of their real-world counter-parts. Each model is defined by a current/voltage equation.
The traditional method for simulation of general electrical circuits and mixed-domain systems is given as a flow chart in FIG. 1. This involves mapping the circuit connectivity information 1001 and the model description 1002 into a set of circuit equations 1003. This is usually done through a formulation, known as modified nodal analysis (MNA). The MNA formulation is well known, and is documented in several references such as in, J. Vlach et. al., “Computer Methods for Circuit Analysis and Design”, New York: Van Nostrand, 1983. These MNA equations are generally in the form of nonlinear differential algebraic equations (NDAEs) and can consist of large number of variables. Such a set of equations is typically represented in the matrix form, where the order of the associated matrices can be very large. For the purpose of illustration, a nonlinear circuit can be represented using the modified nodal analysis in the form of NDAEs, asGx(t)+C{dot over (x)}(t)+f(x(t))=b(t)   (Eq. 1)
where x(t) contains the states of the circuit as a function of time, {dot over (x)}(t) is the time derivative of these states, G is the modified nodal conductance matrix, C is the modified nodal susceptance matrix, f(x(t)) describes the nonlinear devices, while b(t) specifies the voltage and current sources in the circuit.
A time-domain simulation to determine the behavior of a circuit over a time-span of interest involves solution of the above described NDAEs, which typically consists of the following sequence (FIG. 1):
a) Perform a DC analysis (at t=0) 1004 using the initial conditions of the circuit, and obtain the DC solution. The circuit equations for DC analysis generally take the form of nonlinear algebraic equations and are solved using Newton-Raphson iterations.
b) Let tn be the previous time point and set n=0 & tn=0 1005. Determine the next time point for solution, by selecting an appropriate time-step size ‘h’. Let tn+1 be the current time point with tn+1=tn+h 1006.
c) Apply an appropriate integration method (such as Trapezoidal Formula, Backward Euler Formula or any Multi-Step Integration Formula, etc.) 1007 to translate the given set of NDAEs to a set of nonlinear algebraic equations 1008 at the particular time point of interest.
d) Solve the set of nonlinear algebraic equations using Newton-Raphson iterations 1008, until the solution convergence is achieved. This step involves solving the large matrix using LU decomposition and forward/backward substitutions at each iteration. This step becomes a highly time consuming process, when the order of the matrix is large.
e) Check if the predefined truncation error tolerance εte is satisfied 1009. If the error is not satisfied, modify the step size and repeat steps c) to e). If the error is satisfied, check if the maximum time point of interest has been reached 1010. If not, set n=n+1 1011, move to the next time point and repeat steps c) to e), until the time-domain results are computed for the desired time-span of interest.
For the purpose of illustration, if the Backward Euler integration method is applied on the NDAEs represented by (1), then the corresponding difference equation (nonlinear algebraic) is given by
                                                        (                              G                +                                  C                  h                                            )                        ⁢                          x                              n                +                1                                              +                      f            ⁡                          (                              x                                  n                  +                  1                                            )                                      =                                            C              h                        ⁢                          x              n                                +                      b                          n              +              1                                                          (                  Eq          .                                          ⁢          2                )            
where xn=x(tn), xn+1=x(tn+h), bn+1=b(tn+h). To solve this set of nonlinear algebraic equations, algorithms such as the Newton-Raphson method are used, where this equation is reformulated in terms of linearized equations and solved iteratively using LU decomposition and forward/backward substitution (at each iteration). For the purpose of illustration of the steps, let
                              Φ          ⁡                      (                          x                              n                +                1                                      )                          =                                            (                              G                +                                  C                  h                                            )                        ⁢                          x                              n                +                1                                              +                      f            ⁡                          (                              x                                  n                  +                  1                                            )                                -                      b                          n              +              1                                -                                    C              h                        ⁢                          x              n                                                          (                  Eq          .                                          ⁢          3                )            
Let the Jacobian matrix, constructed using the derivative of Φ(xn+1) be represented by
                    J        =                              ∂                          (                              Φ                ⁡                                  (                                      x                                          n                      +                      1                                                        )                                            )                                            ∂                          x                              n                +                1                                                                        (                  Eq          .                                          ⁢          4                )            
Solution at the current time point, xn+1=x(tn+h), involves solving nonlinear algebraic equations of (3) using Newton-Raphson (NR) iterations 1008. The process of NR iterations is detailed in the flow chart of FIG. 2. To begin with, the solution is initialized to the response at the previous time point 2001 during a transient analysis (or to the DC initial conditions during a DC analysis 1004). In each iteration 2002, computing the NR update (Δx) requires forming an updated Jacobian matrix (J) and solving the corresponding equations 2003. For example, at (j+1)th NR iteration, the Jacobian and the right hand side are evaluated based on the solutions from the jth iteration as:
                              J          ⁢                      ❘                          x                              n                +                1                                            (                j                )                                              ⁢                      Δ            ⁢                                                  ⁢            x                          =                              -            Φ                    ⁢                      ❘                          x                              n                +                1                                            (                j                )                                                                        (                  Eq          .                                          ⁢          5                )            
Equation (5) is typically solved using LU decomposition, forward and backward substitution. For example, matrix J can be decomposed into product of two triangular matrices, a lower triangular matrix (L) and an upper triangular matrix (U) (the process is popularly known as LU decomposition), in the form:J=LU  (Eq. 6)
which is a highly time consuming process when the order of the matrix is large. Subsequently, the NR update (Δx) is obtained using forward and backward substitutions as
                                          L            ⁢                                                  ⁢            Z                    =                                    -              Φ                        ⁢                          ❘                              x                                  n                  +                  1                                                  (                  j                  )                                                                    ;                                  ⁢                              U            ⁢                                                  ⁢            Δ            ⁢                                                  ⁢            x                    =          Z                ;                            (                  Eq          .                                          ⁢          7                )            
Once the (Δx) is obtained, the updated solution at the (j+1)th iteration, xn+1(j+1), is obtained 2004. Subsequently, the updated solution is checked for convergence by comparing it against a predefined NR iteration error tolerance, εnr, 2005. If the error is not satisfied, set j=j+1 2006, and move to the next iteration. If the error is satisfied, set the final solution 2007 for the time point tn+1 as xn+1=xn+1(j+1).
In a typical transient analysis involving thousands of time points (note that each time point solution requires several NR iterations), evaluation of the NR update (Δx) using (5) & (6) constitutes the bulk of the computational cost.
An additional factor that influences the computational cost of the solution using (5), (6) and (7) is the sparsity (which corresponds to the number of zero entries in the matrix) of the matrix. In general, computational cost is proportional to Nσ where N represents the original matrix size and σ depends on the sparsity of matrix J. For typical circuits, σ varies from 1.3 to 2.4. Obviously, larger the N, larger is the computational cost. For modern electronic circuits, N can be very large, in the range of millions. Hence it is desirable to approach the circuit simulation problem by dividing the original circuit into several smaller subcircuits (say of size Ns, Ns<<N), and solving each subcircuit independently (to exploit the fact that Nsσ<<Nσ) and combining the subcircuit results to get the solution of the original circuit. This has been attempted through several previous approaches, such as waveform relaxation, as in, J. White et. al., “Relaxation Techniques for the Simulation of VLSI Circuits”, Norwell, Mass.: Kluwer, 1987, and domain-decomposition techniques, such as described in, N. Frohlich et. al., “A new approach for parallel simulation of VLSI Circuits on a transistor level”, IEEE transactions on circuits and systems-I, June 1998, pp. 601-613.
Waveform relaxation techniques approached the large circuit and parallel computing problem by independently computing the entire waveform for subcircuits over a time duration of interest assuming some initial waveforms. These waveforms were iteratively refined to obtain the final waveform for the entire circuit. The main issue with the waveform relaxation approaches is the convergence, which is not guaranteed in the case of general circuits. This has inhibited the use of waveform relaxation in general circuit simulators.
In the case of domain-decomposition based techniques, given the system of original circuit equations (Ax=b), partitions and interface variables are identified. These interface variables are extracted as a subset of the MNA (modified nodal analysis) variables. While ordering the MNA, these interface variables are ordered towards the end. This leads to a domain-decomposed form for the original circuit, with block diagonal matrices for the subcircuits and two border rectangular (one vertical and one horizontal) matrices comprising interface variables as,
                                          [                                                                                A                    1                                                                                                                                                                                                                                                                        F                    1                                                                                                                                                                                                            A                    2                                                                                                                                                                                F                    2                                                                                                                                                                                                                                                                                ⋱                                                  ⋮                                                                                                  M                    1                                                                                        M                    2                                                                    …                                                  C                                                      ]                    ⁡                      [                                                                                x                    1                                                                                                                    x                    2                                                                                                ⋮                                                                                                  x                                          k                      +                      1                                                                                            ]                          =                  [                                                                      b                  1                                                                                                      b                  2                                                                                    ⋮                                                                                      b                                      k                    +                    1                                                                                ]                                    (                  Eq          .                                          ⁢          8                )            
where Aj represents an individual domain (i.e., partition) and Mj, Fj describe the MNA variables that are interfacing Aj to other subcircuits and ‘C’ represents the interconnecting matrix linking all subcircuits. Each partition is solved independently and synchronized using a master process at each Newton iteration. However, the above approach met with little success for parallel computation, mainly due to the fact that the matrices Mj and Fj contain the original MNA variables and are non-binary matrices. Consequently, solution cost of interconnecting equations (solved on a master processor) as well as the communication cost among slaves and the master processor grow rapidly with the increasing number of partitions. This causes poor scalability with the increasing number of processors and partitions. Also, in many cases great care needed to be taken while removing a node from a domain (in order to move it to the border interface matrix), since it could severely affect the solvability of an individual domain, which in turn hampered in the process of efficient partitioning.
Hence, what is needed is a method for accurate and fast transient analysis of large circuits and mixed-domain formulations that effectively partitions the given problem while providing a mechanism requiring minimum computational/communication cost to synchronize the solution among different partitions/processors.